function [ y , t, out ] = simulaceLin( model, linModel, inicializaceMPC, k_max, Ts, x0, u0, param )
%[ y , t, out ] = simulace( model, linModel, inicializaceMPC, k_max, Ts, x0, u0, param )
u=u0;
x(:,1)=x0;
outputSize=[];
t=[];
[A,B,C,D]=linModel(x(:,1),u,param);
system=c2d(ss(A,B,C,D),Ts);
sys.A=system.A;
sys.B=system.B;
sys.C=system.C;
sys.D=system.D;
sys.Ts=system.Ts;
for k=1:k_max;
    %errx=x(:,k);
    %save err errx u k
    %[A,B,C,D]=linModel(x(:,k),u,param);
    if(isempty(outputSize))
        outputSize=size(C,1);
    end
%     if(any(eig(A)>1))
%         x_end=x(:,end);
%         %disp(['vlastni cislo vetsi nez 1: ',num2str(max(eig(A)))]);
%         %save posEig x_end u k
%         %break;
%     end
    mpc=inicializaceMPC(sys,x(:,k),u,k, param);
    [u,up,vecp]=MPCNew(mpc);
    if(isempty(mpc.ref))
        out.ref(:,k)=zeros(outputSize,1);
    else
        out.ref(:,k)=mpc.ref(1:outputSize);
    end
    out.u(:,k)=u;
    out.up(:,:,k)=up;
    out.vecp(:,:,k)=vecp;
    [dState,y_k]=model(x(:,k),u,param);
    [dState2,~]=model(x(:,k)+Ts*dState,u,param);
    x(:,k+1)=x(:,k)+Ts/2*(dState+dState2);
    y(:,k)=y_k;
    t(k)=k*Ts;
end
out.x=x;
end

